

% Hei some nice scripts for auto open, download urls
% you can even zip it later !

fid = fopen('c:\manoj\m\ocean\table.dat','rt');

for i = 1:86,
    [fname,c] = fscanf(fid,'%s',1);
    [id,c] = fscanf(fid,'%d',1);
    dc = num2str(id);
    urlwrite(['http://sio.ecco-group.org/cgi-bin/nph-dods/datasets/glb_1x1_92_02_iter69/iter69t.cdf.ascii?TEMP[' dc '][0:1:22][0:1:159][0:1:359]'],[fname 't']);
end;
%--printing fileout up (not a great deal !)
s = dir('c:\manoj\ocean\60x180\monthly\');
fid = fopen('16x180.bat','wt')
for i = 1:length(s)-2,
fprintf(fid,'tides %s %s\n',s(i+2).name,strrep(s(i+2).name,'.ind','.out'));
end;


%------- divergenec on spherical cordi=nates

sinlat=1./(6378.137*1000)*sin(X(1,:)*pi/180);
sinlat1 = repmat(sinlat,[360,1]);
fac = repmat(sin(X(1,:)*pi/180),[360,1]);
[px,junk] = gradient(U);
[junk,py] = gradient(V.*fac);
DIV = px.*sinlat1+py.*sinlat1;


%---------------------------------------------

%topo silly
load topo;
% 0-359
[c,h] = contour(repmat([0.5:359.5],[180,1]),repmat([-89.5:89.5]',[1,360]),topo,[0 0],'k');
%-180 to 180
[c,h] = contour(repmat([-180:179],[180,1]),repmat([-90:89]',[1,360]),circshift(topo,[0,180]),[0 0],'k');

%--------------------------------------------------------------------------

%projection of one vetor to another (dot product)
% here the anomolus vector field b~n is projected on to the direction of
% main field B

for ij = 1:132,
    bscn(ij,:,:) = dot([bpn(ij,:,:);btn(ij,:,:);brn(ij,:,:)],B(1:3,:,:),1)./B(4,:,:);
end;


%for reading fortran unformatted file wrote to prospect
MyPath = 'C:\manoj\ocean\150x300\monthly\';
nth=150;
nph=300;

Month = ['Jan';'feb';'MAR';'apr';'may';'jun';'JUL';'aUG';'sep';'Oct';'NoV';'deC'];
Month = lower(Month);
Year = ['92';'93';'94';'95';'96';'97';'98';'99';'00';'01';'02'];
ij=1;
nfiles=0;
for i = 1:11,
    for j = 1:12,
        fname = [Month(j,:) Year(i,:) '.out'];
        disp([MyPath fname]);
        fid = fopen([MyPath fname],'r');
        if fid == -1,
            flag=0;
        else,
            flag=1;
            nfiles = nfiles+1;
        end;
        
if flag == 1,
fread(fid,12,'int16');
xx = fread(fid,[150,300],'double');
brn(ij,:,:) = circshift(xx,-1);
xx= fread(fid,[150,300],'double');
bra(ij,:,:) = circshift(xx,-3);
xx = fread(fid,[150,300],'double');
btn(ij,:,:) = circshift(xx,-5);
xx = fread(fid,[150,300],'double');
bta(ij,:,:) = circshift(xx,-7);
xx = fread(fid,[150,300],'double');
bpn(ij,:,:) = circshift(xx,-9);
xx = fread(fid,[150,300],'double');
bpa(ij,:,:) = circshift(xx,-11);
fclose(fid);
end;
ij=ij+1;
    end;
end;

% tha = -(90-acos(1-([1:nth]-0.5)*2/nth)*180/pi); % equal area 
% phi = (-pi+([1:nph]-0.5)*2*pi/nph)*180/pi;

% getting the main field (internal part) B[x,y,z,f]
B = get_pomme(200,tha+90,phi); %fday == 200 just arbitrary - secular varaition is neglig/comp main field

bscn = dot([bpn;btn;brn],B(1:3,:,:),1)./B(4,:,:);

X(1,:,:) = bpn;
X(2,:,:) = btn;
X(3,:,:) = brn;

%20 June 2005 The world geomagnetic observatories table is available in
%c:\geomag\gmagobs.xls. It can be read by xlsread

worldmap([8 15],[74 81]);
load c:\manoj\geomag\figures\out.dat
lat = reshape(out(:,1),[515,648]);
long = reshape(out(:,2),[515,648]);
elev = reshape(out(:,3),[515,648]);
pk = elev;
surfacem(long,lat,pk);
colormap(flipud(gray))
[a,b] = xlsread('c:\manoj\geomag\figures\magnetometer-sites.xls')
A= char(b(:,3))
plotm(a(:,1),a(:,2),'r.','MarkerSize',30)
h=textm(a(:,1),a(:,2)+0.1,A)
get(h)
get(h(1))
set(h,'FontSize',16)
[a,b] = xlsread('c:\manoj\geomag\figures\magnetometer-sites.xls')
A= char(b(:,3))
plotm(a(:,1),a(:,2),'r.','MarkerSize',30)
h=textm(a(:,1),a(:,2)+0.3,A)
get(h)
get(h(1))
set(h,'FontSize',16)
whos
imagesc(pk)

% the following scripts were written after went mad by Alexei's out put
% format This stupid script wades through binary junk to search for useful
% data 22. Aug 2005
flag = 1;
nbytes = 1;
while flag == 1,
    fid = fopen('c:\manoj\ocean\x3dg\1x1_winter_2002_01.res','r','b');
    fread(fid,nbytes,'char');
    nbytes = nbytes + 1;
    d = fread(fid,[180,360],'float');
    plot(mean(d));
    pause;
    fclose(fid);
    fprintf('nbytes = %d\n',nbytes);
end;

% The following scripts are used to read a res file

fid = fopen(['c:\manoj\temp\1x1_winter_2002_new_000.res'],'rb');
fread(fid,518400,'uchar'); %Header includes the gird size, files used for computation etc
[Hpr(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi; %in nT
[Hpi(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
[Htr(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
[Hti(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
[Hrr(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
[Hri(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
[Epr(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*1.E+06; %in mV/km
[Epi(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*1.E+06;
[Etr(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*1.E+06;
[Eti(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*1.E+06;
[Err(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*1.E+06;
[Eri(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*1.E+06;
fclose(fid);

% Reads UPW files

fid = fopen('c:\manoj\projects\cable\apr000.upw','rb');
fread(fid,518400,'uchar'); %Header includes the gird size, files used for computation etc
[Hpr1(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi; %in nT
[Hpi1(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
[Htr1(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
[Hti1(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
[Hrr1(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
[Hri1(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
fclose all;

% jUST TO READ THE RES FILES 21JULY2006 TSUNAMI NUE KOMPUTATION

MyOutPath = 'E:\manoj\projects\cable\result\';
for i = 1:156,
fid = fopen([MyOutPath sprintf('cable%d',i) '.res'],'rb');
fread(fid,518400,'uchar'); %Header includes the gird size, files used for computation etc
[Hpr(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi; %in nT
[Hpi(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
[Htr(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
[Hti(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
[Hrr(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
[Hri(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
[Epr(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*1.E+06; %in mV/km
[Epi(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*1.E+06;
[Etr(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*1.E+06;
[Eti(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*1.E+06;

fclose(fid);
j  = sqrt(-1);
Hr(i,:,:) = Hrr + j*Hri;
Ht(i,:,:) = Htr + j*Hti;
Hp(i,:,:) = Hpr + j*Hpi;
Et(i,:,:) = Etr + j*Eti;
Ep(i,:,:) = Epr + j*Epi;
end;


load c:\manoj\projects\tsunami\20july2006_fields_spec Hr;
%Ht(221:301,:,:) = 0;
Et(241:301,:,:) = 0;

for i = 1:60,
    for j = 1:91,
%         Hr_t(:,i,j) = real(ifft([squeeze(Hr(:,i,j)).' fliplr(squeeze(Hr(2:end-1,i,j))')]));%mult with 2 to remove 1/2 factor in ifft
%         Hp_t(:,i,j) = real(ifft([squeeze(Hp(:,i,j)).' fliplr(squeeze(Hp(2:end-1,i,j))')]))/600;%div with 600 to remove 1/N factor in ifft
%         Ht_t(:,i,j) = real(ifft([squeeze(Ht(:,i,j)).' fliplr(squeeze(Ht(2:end-1,i,j))')]))*2;%mult with 2 to remove 1/2 factor in ifft
          Et_t(:,i,j) = real(ifft([squeeze(Et(:,i,j)).' fliplr(squeeze(Et(2:end-1,i,j))')]))*2;%div with 600 to remove 1/N factor in ifft
%         Ep_t(:,i,j) = real(ifft([squeeze(Ep(:,i,j)).' fliplr(squeeze(Ep(2:end-1,i,j))')]))/600;%div with 600 to remove 1/N factor in ifft
%Note in the earlier computation, since the spectra was scaled by N, the
%fields were multiplied by 1200
    end;
end;


% Following scripts used to generate the dip_latitude and eej-factro 10 Sep
%2006
load c:\manoj\geomag\dip_eq.mat lat111 long111;
load c:\manoj\projects\eej\OBS_array.mat OBS;
%[a,b] = xlsread('c:\manoj\projects\eej\ETT-HYB-CHAMP.xls');
data = load('c:\manoj\projects\eej-induction\GROUND_5.TXT');
a(:,1) = data(:,1);
a(:,2) = data(:,2)+data(:,3);
[aa,bb,cc] = unique(long111);
obs_pairs = [2,1;2,3;2,13;2,16;2,17;2,18;4,1;4,2;4,3;4,13;4,16;4,17;4,18;5,6;7,8;9,10;9,14;15,10;15,14;...
18,1;18,3;18,13;18,16;18,17;13,1;13,3;19,20];
for i = 1:20, if OBS(i).long > 180, OBS(i).long = -(360-OBS(i).long); end;
OBS(i).dip_lat = OBS(i).lat-interp1(long111(bb),lat111(bb),OBS(i).long);
OBS(i).eej_factor_new = interp1(a(:,1),a(:,2),0)./interp1(a(:,1),a(:,2),OBS(i).dip_lat);
end;

obs_pairs = [2,3;4,1;5,6;19,20;15,14];
fid = fopen('c:\manoj\temp\eej_slope.txt','wt');
for i = 1:length(obs_pairs),
fprintf(fid,'%s-%s ', OBS(obs_pairs(i,1)).code,OBS(obs_pairs(i,2)).code);
fprintf(fid,'%5.2f %5.2f %10.5f %10.5f\n',OBS(obs_pairs(i,1)).dip_lat,OBS(obs_pairs(i,2)).dip_lat-OBS(obs_pairs(i,1)).dip_lat,...
    OBS(obs_pairs(i,1)).eej_factor,OBS(obs_pairs(i,1)).eej_factor_new);
end;


%To plot the eej-model-ground.fig
[p,s] = polyfit(data(:,1),data(:,2)+data(:,3),20);
hold on;
for i = 1:length(obs_pairs),
      yy = polyval(p,OBS(obs_pairs(i,1)).dip_lat);
      plot(OBS(obs_pairs(i,1)).dip_lat,yy,'ks');
      text(OBS(obs_pairs(i,1)).dip_lat,yy,OBS(obs_pairs(i,1)).code);
end;
    
% Note with the dip_latitude as determined from EEJ loc form stefan's new
% file, some discripancies have been noted with refernce to an earlier
% similar estimates (given in slop_dip.xls


load coast;
data = flipud(Etr);
worldmap('world');
h=contourfm(data,[1,89.9,0.5],20);
geoshow(lat, long,'color','k');
caxis([-.05 .05]);
title('Im B_{\theta}');
colorbar;
